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Abstract. We argue that when the number of spins N in the SK model is finite, 
the Parisi scheme can be terminated after K replica-symmetry breaking steps, where 
K{N) oc N^/^. We have checked this idea by Monte Carlo simulations: we expect 
the typical number of peaks and features R in the (non-bond averaged) Parisi 
i-rt ■ overlap function Pj{q) to be of order 2K{N), and our counting (for samples of size 

^H ' A^ up to 4096 spins) gives results which arc consistent with our arguments. We 

can estimate the leading finite size correction for any thermodynamic quantity by 
finding its K dependence in the Parisi scheme and then replacing K by K{N). Our 
predictions of how the Edwards- Anderson order parameter and the internal energy of 
^ . the system approach their thermodynamic limit compare well with the results of our 

lO ! Monte Carlo simulations. The A^-dependence of the sample-to-sample fluctuations of 

'nI ■ thermodynamic quantities can also be obtained; the total internal energy should have 

sample-to-sample fluctuations of order N^^^, which is again consistent with the results 
of our numerical simulations. 
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PACS numbers: 75.50.Lk, 75.10.Nr, 75.40. Gb 

1. Introduction 

The Sherrington-Kirkpatrick (SK)[T] model of spin glasses has been the subject of 
hundreds of papers. It is the model for which mean-field theory becomes exact in 
the thermodynamic limit (i.e. when A^, the number of spins in the model, becomes 
infinite). Parisi's replica symmetry breaking (RSB) solution[2] is now known to be the 
correct mean-field solution|3]. Extensive studies, mostly numerical, have been made 
of the model at finite A^ values. Analytically the determination of the properties 
of the model at finite A^ - the finite-size corrections - is a much more challenging 
task in the low-temperature phase than finding the mean-field theory. At finite A^ 
all the loop corrections to the mean-field solution need to be considered. Because 
of the massless modes present in the low-temperature phase [4] each term in the loop 
expansion is infinite. Hence a direct perturbative approach is impossible. A similar 
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situation applies in finite dimensional spin glasses, already for the bulk term, when 
the dimension d is smaller than six[5], so that one could hope at least that experience 
gained in studying finite size effects in the SK model might be relevant to spin glasses 
in physical dimensions. 

Unfortunately we have been unable to find any systematic theoretical treatment 
of the finite-size problem. However, we have managed to obtain insights into it by 
examining the structure of the Parisi overlap probability distribution function Pj{q) (i.e. 
the non-averaged overlap probability distribution function) at finite A^ values. Pj{q) is 
defined as the probability that the overlap of the spins in two copies of the system with 
the same realization of the quenched disorder Jij is equal to q, i.e. 



P,(g)^ %--Va,r,) , (1) 




where the Hamiltonian of the two-copy system is 

'H = - J2 JiMi^^i + ^i^j) ' (2) 

and the sum runs over all the pairs ij of sites in the system. The thermal average (■ ■ ■) in 
Eq. ([1]) is taken over the Boltzmann weight associated with all the possible values (±1) 
of the Ising spins ai and Tj. It has been known for many years that the function Pj{q) 
is very different for different realizations of the bonds. In particular it contains a very 
variable number R of peaks, humps or shoulders. In this paper we shall systematically 
study the distribution of R and the dependence of its average on the number of spins 
A^. We shall give numerical and analytic arguments that the mean number of features 
R increases as N'^, with /i = 1/6 and that 6R, the width of the distribution of R, is 
N independent for large N. The next step of our approach is to argue that the Parisi 
replica symmetry breaking scheme, which involves K levels of symmetry breaking (where 
in order to achieve a stable solution K has to be taken infinite in the thermodynamic 
limit), is stabilized at finite A^ at a value K{N) by self-energy contributions (whose A^ 



dependence is estimated in Appendix Appendix A). As a consequence we can estimate 



K for a given system size and because R = 2K (see Sec. [8]), we can understand the size 
dependence of the number of peaks/features in Pj{q)- 

The next step towards predicting the exponents which give the leading A^ 
dependence of the corrections to the thermodynamic limit of quantities such as the 
internal energy per spin e is simply to use the RSB scheme to compute the dependence 
of the quantity on K. For example e = ep-\- 0{K~^), where cp denotes the value of the 
internal energy in the infinite K limit [6]. Our prescription for evaluating the exponents 
of the leading finite size corrections is to set K = K{N) ~ N^^^; this implies that the 
leading finite size correction to the thermodynamic limit of the internal energy per spin 
should be of order N~'^^^. Since arguments of this type do not have the strength of 
a theorem and can only be suggestive of the possible behaviour, we have checked our 
arguments with extensive Monte Carlo simulations. We have computed many quantities 
at a number of different values of the temperature. Our results for the internal energy 
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are reported in Sec. O That data strongly supports the value of 2/3 predicted by our 
approach for the exponent of the leading finite-size correction. 

Similarly, the Edwards-Anderson order parameter qea at finite K differs from its 
infinite K form by a term of order K~^. Hence we predict that the finite size shift of 
qea should be of 0{1/N^^^), and we present numerical evidence for this behaviour in 
Sec. El 

Our approach can be used to investigate the sample-to-sample fluctuations of any 
quantity by relating them to the sample-to-sample variation in the number of features 
in Pj{q)i SR. For the internal energy we shall find in Sec. [7] numerical evidence 
consistent with this approach, together with a discussion of the behaviour of the sample- 
to-sample fluctuations in the critical regime and in the high-temperature phase. Our 
basic prediction is that the sample-to-sample fluctuations in the total free energy of a 
system of A^ spins are of order A^^ where the exponent T = // = 1/6. There have been 
numerous attempts to determine this exponent, both numerically and analytically, and 
we review them also in SecJTl 

Because the peaks/features in Pj{q) are caused by the overlap of pure states, in 
particular those states whose free energies are of order ksT from that of the lowest free 
energy state, one can relate the number of these pure states to the number of peaks 
R using the relation R = 2K. This connection is simplified because of the ultrametric 
organization of states in the SK model and the details of the argument are given in 
Sec. [HI In Sec. |9] we discuss the relation of these ideas with the behaviour of finite 
dimensional spin glasses. 

2. Theoretical framework 

Our Monte Carlo studies of the Parisi overlap probability distribution function Pj{q) 
for systems of A^ spins (with A^ up to 4096) show that the number R of peaks/features 
is usually quite small, and that it increases only slowly with A^, apparently as i? ~ A^'^, 
with fi ^ 1/6. Our approach to the study of finite size effects in the SK model is to 
argue that R, the average number of such peaks/features for a system of size A^, can be 
connected to a truncation of Parisi's RSB scheme at its i^th step, with R = 2K{N). 

The Parisi scheme at the i^-th level of RSB parametrizes the bond-average of Pj{q), 
P{q), by a series of delta functions at various values of q, viz gi, g2, • • • , qK] 

P{q)=Y.a,5{q-q,). (3) 

The weights of the delta functions Oj and their positions qi are the variational parameters 
that one optimizes to obtain the Parisi solution. In the thermodynamic limit, where 
A^ goes to infinity, a Parisi RSB solution with K > 1 is only stable if K is taken to 
infinity. We argue that in finite size systems the self-energy corrections to the Parisi 
solution can stabilize an RSB solution with a finite value of K, and we will argue 
that R = 2K ~ A^^. (In zero field Pj{q) = Pj{—q) so the number of peaks/features 
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R = 2K. However, if gi just happens to be zero, i.e. there is a peak at the origin, then 

R = 2{K-l) + l = 2K -1). 

Consider the single-valley replicon correlation function G{i{i,j) = {SiSj)^. At 
wavevector k its Fourier transform takes the form described in Ref . [1] , and at Gaussian 
order, Gjiik) = 1/k'^, both for T < Tc and T = T^. (Strictly speaking in the SK model 
the only possible value which k can take is zero, but we will find it useful to consider 
non-zero values of k). Right at A; = 0, Gr{0) is infinite in the thermodynamic limit. 
For finite N, the self-energy corrections neglected at Gaussian order will be shown in 
Appendix Appendix A to produce a divergence growing as N'^^'^. Schematically 



SO that the self-energy S/j is of order 1/A^^/^. 

Now for finite values of K in the Parisi RSB scheme, the Gaussian propagator is 
unstable and behaves as [6] 

GR{k) = / ,. , (5) 

^ 3 {2K+1)'^ 

in the regime near the transition temperature T^ where t = 1 — T/T^ is small. The 
instability at /c = only disappears when one takes the infinite K limit. Our basic 
idea is that for finite A^ this instability can be removed by the stabilizing effect of the 
self-energy S/j. Then if E/? = c/N^^^ stability will be achieved when 

In other words, when K = K{N) ~ tN^'^, there will be no need to break the 
symmetry further (at least to achieve stability). This would explain why the number of 
peaks/features in Pj{q) increases as A^^'^ (see Fig. [1]). 

Our procedure to determine the finite size corrections to scaling of any 
thermodynamic quantity proceeds in a similar fashion. First one obtains from RSB 
calculations the K^^ approximation for the quantity. Thus the free energy per spin 
below but near T^ is to order t^, and at large values of /C [6] 

■' V6 24 120 / 360 \kJ ^ ' 

To estimate the A^ dependence of the finite size corrections we replace K by tN^^^. 
This gives a term in A/ which scales as t/N'^^^, which is in excellent agreement with 
numerical studies [7]. Just as the self-energy corrections to Eq. (^ change the sign of 
Gr{0), we would expect that the higher loop corrections to the free energy will also 
change the sign of this correction, but not its A^ dependence. 

A similar argument can be given for other quantities. The additional terms in the 
internal energy per spin below T^ at order K in the RSB procedure arejS] 
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Figure 1. Scaling plot oi E{R) (the average number of peaks/features determined by 
visual inspection of the individual Pj{q)) as a function of N, for T — 0.4. The curve 
is the best fit to the form E{R) = a + bN" with c = 0.17 ± 0.14. 



Substituting as before tN^^^ for K, the finite size corrections to the internal energy 

would be expected to be of order 1/N'^/^. 

The Edwards- Anderson order parameter [6] is to order K 

2 
qj^j^ = t + t^ -t^ , (9) 

correct to order t^. It is thus to be expected on substituting for K that the finite size 
corrections to qea are of order 1/A^^/^. If one defines qea for finite-size systems as the 
value of q at which the Parisi overlap function P{q) peaks, then such an N dependence 
is in excellent agreement with both existing numerical and theoretical arguments[S]. We 
postpone to Section [6] the comparison with the results of our numerical analysis of the 
scaling behaviour of qEA- 

Our approach can be extended to determine the N dependence of sample-to-sample 
fluctuations of, say, the internal energy or the free energy. In Ref. [9] it was shown that 
the variance of the sample-to-sample fluctuations of the extensive free energy 5F varies 
as a quantity — J(0), which at Gaussian order has the property —J{k) ~ fc"^. At flnite 
RSB of order K, an exact expression for this quantity was given in Ref. [lO] . As shown in 
Appendix Appendix B it can be evaluated when k'^ (which is originally a wave vector) 
is replaced by the self-energy. One gets that 

- J ^ N'I'fit) , (10) 



where f{t) is some known function (see Appendix Appendix B). This shows that the 
variance of 5F is of order N^/^ ^ with typical fluctuations being of order N^^^ . 

From this result one can compute the sample-to-sample fluctuations of R. To do 
this we shall suppose that the sample-to-sample fluctuation of K is of order 5K. Then 
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Eq. ([7]) implies that the sainple-to-sample fluctuation of the extensive free energy SF 
behaves as 

6F ~ Nt^ (^) 5K . (11) 

Given that 5F ^ N^/^ , it follows that 5K is of 0(1), that is, independent of A^. 

We have determined the sample-to-sample fluctuations of the internal energy in 
the course of our simulations, and their A^ dependence can be predicted by extension 
of these arguments. From Eq. ([8]), the sample-to-sample variation of the full internal 
energy 5U is 

5U ^ Nt^ (^\ 5K . (12) 

Substituting for K and 5K^ it follows that 5U ~ t~^N^^^. These sample-to-sample 
fluctuations appear to diverge at T = Tc, but Eq. [12] only holds in the RSB region, 
which is outside the critical regime, (which is defined by the limits A^ -^ oo, t ^ 0, with 
iVt^ fixed HH). 

All these arguments are intuitive rather than rigorous. As a consequence we have 
attempted to check them by numerical simulations of the finite size SK model. 

3. The Monte Carlo simulation 

We have based our analysis on a large set of numerical data produced by the large scale 
parallel tempering simulation of Ref. [12] and [13], supplemented by a new large scale 
simulation for lattices with A^ = 2048 spins. 

The quenched random couplings of our system can take the two values ±1 with 
equal probability; the use of such binary couplings allows to write computer codes that 
run much faster than, say, when using quenched random couplings assigned under a 
Gaussian distribution. We assume that the interesting leading scaling behaviour is the 
same, for example, when using binary or Gaussian couplings. 

We report in Table [1] the relevant parameters of our numerical simulations. The 
temperatures allowed to the parallel tempering steps are in the range T G [0.4, 1.3]. 

A parallel tempering sweep consists of one Metropolis sweep (all spins are updated 
in lexicographic order) followed by a temperature exchange sweep (we try to exchange 
adjacent values of T in sequential order). The balance between the number of sweeps 
performed for each disorder sample and the number of disorder samples included 
has been chosen cautiously in order to avoid any possible bias due to a non-perfect 
thermalization; we have chosen a safe compromise favouring, at fixed amount of 
computer time, the number of sweeps over the number of disorder samples. We have 
checked the quality of thermalization by monitoring for example the value of q^ as a 
function of the Monte Carlo time, starting from an ordered initial spin configurations 
(all spins equal to one). For all values of T and A^ the disorder averaged data do not drift 
appreciably already after a couple of thousands sweeps, i.e. far before we start taking 
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Table 1. The relevant parameters of our numerical runs: number of sites N , number of 
parallel tempering sweeps used for measurements Nmeas , number of parallel tempering 
sweeps used for thermalization Nequi, and number of disorder samples Nj. 



measurements. A second important test is provided by the symmetry of the individual 
Pj{qys that is very good for most samples (See Sect. H]). 

In the rest of this note we will denote by E{- ■ ■) the average over the quenched 
disorder, U = Ne the total internal energy and and 6U = NA its standard deviation, 



with N^A^ = E{U^) - EiUf = N^(^E{e^) - E{e) 



4. The structure of Pj{q) 

As mentioned before, the function Pj{q) is very different for different realizations of 
the bonds. Fig. [2] shows eight such distributions for the lowest temperature value 
(T = 0.4) of the largest system (A^ = 4096) we have simulated. The symmetry of 
the plots under inversion of the overlap is excellent: one can see that even very small 
peaks appear with their reflected counterpart, and this is a remarkable check of good 
thermalization. The only mild asymmetries one can see concern the peaks heights, 
connected to the population of the different "pure states to be" , that is a very difficult 
quantity to estimate by Monte Carlo integration (just think about the two peaks in the 
magnetization distribution for the usual Ising model in three dimension below the Curie 
temperature). Each of the Pj exhibits a given number Rj of features (well formed peaks, 
humps, shoulders on the side of a peak, and so on) that we are interested to determine 
in order to compute its disorder expectation value E{R) and the scaling behavior of 
E{R) with A^. 

Peaks/feature counting is not at all a trivial issue when dealing with noisy data. 
The first assumption must be that the statistical accuracy of the data set is good enough 
not to hide important features (in this respect it is possible that further improvements 
to the Monte Carlo scheme could help: one could check for example if using multi- 
overlap algorithms [m [15] could be of help). Under this assumption we have developed 
a computer code that counts the number of peaks. The first step is based on smoothing 
up the (symmetrized) data for Pj{q)- The second step determines the peaks of the 
soothed data: the peaks are defined as local maxima of Pj{q) where a valley at least 
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Figure 2. Pj{q) for eight different disorder realizations: here N = 4096 and 
T — 0.4000. The symmetry of the plots around g = is a good test of thermalization. 
The number of peaks R quoted above each figure is the value computed by our computer 
program. 



lower than a percentage p of the peak height follows the peak on both sides (we have 
selected p = 90% and we have relaxed the depth condition for valleys that include one 
of the two frontiers of the support, i.e. q = ±1). As a third and last step we impose 
a cutoff on the putative peak heights: we discard any peak of height lower than 1% of 
the highest peak present in the Pj{q) for the given disorder sample. In what follows, we 
will call "automatic peak counting" this procedure to determine features. 

We did not push the coding to include in our automatic peak counting the more 
complex structures which visual inspection spots: an automatic approach to such a 
complex task needs great care to avoid arbitrary choices that could lead to misleading 
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conclusions. 

For example it looks clear that the first plot of Fig. [2] (the upper left figure) is 
characterized by four features, two for g > and the two symmetric ones for g < 0. 
The first feature is the clear peak that also the computer code finds, while the second 
very clear feature is the shoulder on the peak: this shoulder is naturally interpreted as a 
second unresolved peak, too wide to be an isolated feature, but whose presence is very 
clear to the observer. In other words it is clear that a correct analysis of this feature 
would lead to R = A and that the conclusion R = 2 reached by our computer code is 
not careful enough (as an optimistic remark let us add that it is possible that the lazy 
procedure could lead asymptotically to the same scaling behaviour of a more careful 
counting). To be on the safe side, we have carried out the tedious task of looking at the 
first 192 Pj(g)'s for T = 0.4 and estimating by eye the number of features R for every 
graph. This procedure will be called "visual inspection" in what follows. 

We show the behaviour of E{R) as a function of A^ for T = 0.4 in Fig. [1] for 
the visual inspection and in Fig. [3] for the automatic peak counting. The statistical 
errors are estimated from the fluctuations between the disorder samples. Our visual 
inspection gives on average between one and two extra features that are not found 
by the automatic computer counting. The best fits to the form E{R) = a + hN"^ for 
the two cases give the exponents Cyisuai = 0.17 ± 0.14 for the visual inspection and 
Cautomatic = 0.28 ± 0.04 for the automatic count. The difference between the estimated 
error bars mainly reflects the number of disorder samples considered in the two cases. 
Our prediction is c = 1/6 ~ 0.17. Our result for casual is on the top of it but it has a 
huge statistical error, while our result for Cautomatic is not consistent with 1/6 at more 
than two standard deviations (but is of uncertain relevance). It is important to notice 
that here we are dealing with a quantity that grows very slowly and is very small even 
for our largest systems, and that pre-asymptotic effects could be large. 

In Fig. IHwe show 

6R = [E{R^) - E{RfY/^ , (13) 

i.e. the fluctuations of R (as determined by visual inspection), together with the best 
fit to the form 5R = a + hN'^. The best fit is obtained for c = 0.07 ± 0.13: our 
theoretical estimate in Sec. [2] of the value of the exponent c was zero; our numerical 
work is consistent with this estimate, although greater precision is really needed before 
the result can be regarded as definitive. 

In Fig. [5] we give the empirical distribution of the number of features R obtained 
by visual inspection for A^ = 4096 and T = 0.4: the distribution is very wide. It is also 
clear from the figure that even values of R are more common than odd values: this is 
expected, since odd values are only obtained in cases where Pj{q) has a peak in g = 0. 

5. Finite size shift in the energy 

In Fig. [6] we show the internal energy per spin as a function of A^~^/^ for our lowest 
temperature T = 0.4 (that we believe is low enough to be free of effects from the critical 
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Figure 3. Scaling plot of E{R) determined by automatic peak counting, as a function 
of iV, for T = 0.4. The curve is the fit to the form E{R) = a + bN" with c = 0.28±0.04. 
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Figure 4. Scaling plot of 5R determined by visual inspection, as a function of N for 
T ~ 0.4. The curve drawn is a fit to the form 5R = a + bN'^. The best fit is obtained 
for c = 0.07 ±0.13. 
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Figure 5. Distribution of the number of features R obtained by visual inspection for 
N = 4096 and T = 0.4. 
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Figure 6. The internal energy as a function of N ^/■^ for T ~ 0.4. The line is a linear 
fit (using data for N > 256) as a function of N^"^/'^ to the data. 



point). The statistical errors are again estimated from the fluctuations between the 
disorder samples. 

In the same flgure we show a flt to the form cat = Coo + ^ A^~^'^, using the value 
Coo = -0.735110726 from Ref. [E]. The best flt is obtained for A = 0.77 ± 0.01, with a 
X^ of 12 for 4 degrees of freedom. The presence of slow decaying sub-leading corrections 
(the dominant sub-leading contribution behaves like 1/A^, barely faster than A^~^/'^) 
explains presumably why the x^ is larger than the number of degrees of freedom. The 
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Figure 7. N'^/^{eN - eoo) as a function of N ^/^ for T = 0.4. The line is a linear fit 
(all data points are included) to the data, of the form A + B N'^^^. 



importance of the corrections to the leading behavior (together with some statistical 
oddness) shows up in Fig. [3, where we plot N'^^^{eN — Coo) as a function of 1/A^^/^: this 
is a way to focus on the deviations from the leading behaviour. We do believe that the 
three leftmost odd looking data points in Fig. [7] are due to a statistical fluctuation. The 
curve is a fit to N'^^^{eN — eoo) oc 1/A^^/^, which is the form we would expect from a 
1/A^ correction to the internal energy per spin. Based on Figs. El Hand similar plots at 
different temperature values, we conclude that our numerical data are consistent with 
an exponent 2/3 in the whole spin glass phase. We disagree with the conclusions of 
Ref. [T7j that are based on data with 36 < A^ < 196, namely a region that is discarded 
altogether in our fits (see Ref. [18] for a detailed comparison). 



6. The Edwards-Anderson order parameter qea 

We have studied the finite size behaviour of the Edwards-Anderson order parameter 
qea, improving the analysis of Ref. [S]. We define qea on a finite system as the location 
of the maximum of the disorder averaged P{q) = E{Pj{q)). The exact procedure is 
the following: we first symmetrize our data for P{q), then we determine the maximum 
value reached by the function (we call it Pmax) and finally we compute qea by means of 
a quadratic fit of the data in the range of positive q values such that P{q) > 0.95Pmax, 
using the same weights for all data points. This gives us an estimate of qEA^N) that 
is not forced to take discrete values: statistical errors are obtained through a jackknife 
analysis (obviously a jackknife approach would not make sense if qEA was constrained 
to take discrete values). 

We have compared the values we have obtained for qEA to values obtained with 
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Figure 8. The Edwards- Anderson order parameter qEA{N), as defined in tiie text; 
here T = 0.5. The line is for the best fit to the data as a linear function of A^"^/^ (for 
N > 256). 



shorter numerical simulations, and there is an excellent agreement: this strongly suggests 
that the procedure we have used to determine qEA{N) has no appreciable statistical bias. 
We show in Fig. [8]our data for qEA{N) as a function of A^~^/^ for T = 0.5, together with 
our best fit (using data with A^ > 256) that uses the infinite volume result qea = 0.6395 
from Ref. [16]. The value of x^ is 1-3 for 4 degrees of freedom, suggesting that in this 
case sub- leading corrections are not very relevant. The small A^ data exhibit larger 
corrections from the asymptotic behaviour than the ones for the internal energy. This 
is not unexpected since the definition of qea on a finite system is involved: for example 
the intrinsic resolution of the determination of qea in a finite system is 2/N (that is 
close to 0.03 for A^ = 64), a value that is exactly on the scale of the deviations that we 
observe. 



7. The sample-to-sample fluctuations of the internal energy 

We have also analyzed the fluctuations of the internal energy between different bond 
realizations using as a measure 

6U' = N\E{{e)])-E{{e)jr). (14) 

The issue of the scaling behaviour of the free energy fluctuations in the SK model has 
been investigated intensively over the years. Let us define an exponent T by the equation 



6F' = N'{E{{f)])-E{{f)jr)^N'^ 



(15) 



Our theoretical approach suggests that T should be /i = 1/6. The first investigation 
we know of is the numerical work of Ref. [19] at T = which gave (on very small 
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samples) the result T = 0.222 , compatible with T = 1/4. The later theoretical analysis 
of Ref. [20] gives T = 1/6. However, there is a caveat to this conclusion. What is 
effectively calculated is the tail of the probability distribution of the free energy on 
the low-free energy side [21]; in order to get T and one has to assume that the N 
dependence of the fluctuations in the tail equals the one of the standard deviation of 
the free energy. Ultimately these analyses relate the value of T to the order of the 
first non-linear term in the expansion of the replicated free energy in powers of n (the 
number of replicas), which is the n^ term [1]. More recently, several authors, using 
exact or heuristic ground states determination algorithms, have found zero temperature 
values compatible [221 ESI El ESI ESI [27] with T = 1/4, some excluding [22l ESI ESI EE] 
the T = 1/6 value, some not [211 [27]. Analytical arguments in support of the T = 1/4 
value can be found in Refs. [28l |24]. Finally a recent numerical simulation [29] using 
an innovative method obtains T = 2/5 in the low T phase. Clearly, the situation is far 
from being settled. It should be clear that in the numerical approach it is extremely 
hard to distinguish with confidence exponents as close as 1/4 and 1/6 when the range 
of variations of A^ is small (typically one decade), the more so as the functional form of 
the next correction is unknown. A further difficulty is that exact algorithms are limited 
to very small systems and heuristic algorithms are heuristic. 

The above results are for the free energy. The problem for the internal energy 
at finite temperature has never, to our knowledge, been studied numerically; we will 
consider here finite, low temperature values, and this will allow us to analyze both the 
low T phase and the T — > limit. 

Let ej be the energy measured at some temperature T during the numerical 
simulation of a system with a given realization of the disorder J (this is what we call a 
sample of our system), ej comes from the average of N^eas values, and we can use the 
quantity 



(16) 



A\T) = E{{e')^)-E{{e)jy. (17) 



/^J J J 

where Nj is the number of samples, to estimate 

At leading order A^(T) and SU"^ are related through 

^, «^ 2Tg;r^£^ (18) 

where Ti^t is the integrated autocorrelation time for the energy at temperature T, Cn(T) 
is the specific heat and N^eas is the number of parallel tempering sweeps performed 
during the measurement phase of the simulation (we measure the energy after every 
Metropolis sweep of the system). 

It turns out that in our numerical data the second term in Eq. [18] is negligible, as 
we have checked by comparing the estimates of A^(T) in different numerical simulations 
for the same set of disorder couplings (with A^ = 64, 256 and 1024) using the first 200i^ 
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Figure 9. The energy fluctuation SU /N ~ iVA^(T) as a function of T, for different 
values of A^. 



parallel tempering (PT) sweeps and the second lOOOJir PT sweeps after thermalization. 
To the best of our knowledge the autocorrelation time Tj^^ of the parallel tempering 
algorithm has never been measured for the SK model, and our results show that it is 
very small. Our results for A^(T) as a function of T can be found in Fig. [9l 

Fig. [To] shows the scaling behaviour of 6U for our lowest temperature. The leading 
exponent is compatible with the value 1/6 but our statistical (and systematic) accuracy 
is not good enough to allow us to rule out the value 1/4: the main culprits are the data 
points for large systems (mostly N = 4096 and A^ = 2048), and we would need a much 
larger number of samples to have a precise determination of this exponent from Fig. [TUl 
only. 

In the critical region SF ~ f{tN^'^), where SF are the sample to sample fluctuations 
of the total free energy [20] and t = 1 — T/Tc. The sample to sample internal energy 
fluctuations are related to the t derivative of 6F and they scale as N^^^f{tN^^'^), where 
/(x) goes like 1/x at large negative x, see Fig. [11], and is a constant at x = 0, i.e. T = Tc. 
The scaling is excellent in the paramagnetic phase (on the left) and in the spin glass 
phase where tN^^^ is small, namely before the oo-RSB effects start to be important, 
leading to a different behaviour. Fig. [11] shows clearly the two scaling regimes. The 
oo-RSB effects are only present for T < T^ in the regime where Nt^ is large [30] so 
multiple pure states can exist. In the finite-size critical regime one has Nt^ fixed with t 
going to zero. This makes Nt^ go to zero, so that in the critical regime RSB effects are 
absent. 

In Fig. [12], we visualize the scaling behaviour of 6U in a different way; we show the 
exponent 1 — C obtained from a fit of x = SU/{E{e)) to the form oc A^^^'' as a function 
of T. This plot is consistent with the guess that the exponent at T = takes a value of 
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Figure 11. Scaling plot of SU/N^/^ as a function of iA^-^^^. 



1/6: in order to get a value of 1/4 we should have a very complex T dependence of the 
effective exponent. The situation at T < Tc and exactly at T = T^ is more complicated: 
it is possible to see that finite size effects bring down the value of the exponent with 
increasing lattice size, and a scenario where the exponent is 1/6 for all T < Tc (but with 
large finite size corrections) is plausible and consistent with the data. 
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8. The number of pure states 

It is useful to consider how the peaks in Pj{q) arise. Suppose we have just two states, 
state a and its spin reverse A. Then using the definition of Pj{q) in Eq. [H if copy 1 is 
in a and copy 2 is also in a there will be a peak at at +qEA- If both copies are in A, 
there will also be a peak at +qEA- However, if copy 1 is in state a and copy 2 is in state 
A, that will produce a peak at —qEA- the same will be true when copy 1 is in state A 
and copy 2 is in state a. 

Suppose now we have 4 states, a. A, b, B. Because (in the infinite volume limit) 
qEA is the same for all states the overlaps aa, AA, bb, BB are all qEA-, and aA, bB, Aa, 
Bb are all —qEA- The overlaps ab=qi2=AB, and aB=-qi2=Ab, giving 4 peaks in total. 
(If qi2 = we have 3 peaks only). However, the peak at ±gi2 will not in general have 
the same weight as that at qEA- 

Then, from the above, 2 states give 2 peaks and 4 states give 4 peaks (if all involved 
overlaps are large than zero). 

With three states 1, 2 and 3 the effects of the ultrametric organization of states 
start to play a role. Besides the peak at qEA, there could be overlaps gi2, 923 and gi3 
making at most 4 peaks (if all the overlaps are nonzero.). But ultrametricity says that 
either all three q's are equal or 2 are equal, making at most 3 peaks (6 peaks when one 
includes the time-reversed states). 

With four states 1, 2, 3 and 4, besides the peak at qEA, there are the overlaps qu, 
qi3, qu, q23, q2A, qzi, but here ultrametricity limits one to 3 distinct possibilities, making 
4 peaks in total (and 8 when one includes the time-reversed states). 
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There is clearly a pattern here; if states are organized ultrametrically their number 
is equal to the number of peaks of Pj{q) (plus one if there is a peak at g = 0). We 
are thus predicting that the number of pure states grows with A^ as ~ N^^^. These 
pure states are those whose total free energy is of order ksT from that of the lowest 
free energy state. It is only the overlaps of these low-lying states which can produce 
detectable features in Pj{q)- 

Another way of determining the number of pure states at level K is to observe that 
if one starts from the bottom (the leaves) of a genealogical tree, and moves towards the 
ancestors, points that are going up can only meet at the bifurcations of the tree. Hence 
the number of overlaps is equal to the number of levels K in the tree. 

9. Application to a finite number of dimensions 

The exponent T which determines the sample-to-sample fluctuations SF of the free 
energy of the SK model has a significance beyond this model as it appears in the theory 
of the interface free energy of finite dimensional spin glasses. 

In Ref. [28] it was shown by going to one-loop order about the Parisi RSB 
mean-field solution that the variance of the interface free energy associated with the 
change in the free energy on going from periodic to antiperiodic boundary conditions, 
SPp.AP = Fp — Fap, is of the form 



5FlAP = L^f{L/M) + 5F^. (19) 

Here the system is of length L in the z direction, and it is periodic and of length M in 
the transverse d — 1 dimensions. The change from periodic to anti-periodic boundary 
conditions is done by flipping the sign of the bonds in a hyperplane perpendicular to 
the 2;-axis. It follows that 6Fp ap = 0. 6F'^ is the bond-averaged variance of the free 
energy of the SK model containing A^ = LAF^'^ spins. Eq. ( TT9l) is valid at least to one 
loop order, and its form is probably unchanged whenever the loop expansion is possible; 
the loop expansion is well-defined in the low-temperature spin glass state when d > 6, 
but its existence is problematic for rf < 6 (see Ref. [1]). 

The first term in Eq. (fT9|) is of the standard aspect-ratio scaling form [311 [32] . 
where the zero-temperature scaling exponent 6 is equal to 1. When L is of order M, 
SF"^ is of order A^^^, i.e. of order ~ L'^^^ if T = 1/6. This term is not of the standard 
aspect-ratio scaling form; it depends instead on the total number of spins in the system. 
This reflects the fact that in RSB situations domain walls have a fractal dimension dg 
equal to d, i.e. they are space filling. For example in dimension d = 2, where we know 
that we do not have a broken symmetry phase, domain walls are fractal with dg < d. 
The variance of the interface free energy is dominated by the SK like term for all d > 6 
provided that T = 1/6. If T = 1/6, one would expect that numerical studies of the 
defect energy in six dimensions would suggest a value of 6 close to unity; exactly ind = 6 
Boettcher|25] found that ^^l.liO.l. It is possible that when d < Q the standard 
aspect-ratio scaling form will dominate. (Of course, when d < 6, the one-loop expression 
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for the interface free energy will no longer be adequate). Thus the dominant term in 
the L dependence of the interface free energy could have very different forms above and 
below six dimensions; note that d = 6 plays a special role only if the exponent T is 
exactly 1/6. 

For d < 6 the loop expansion about the Parisi RSB state becomes problematical 
but it is clearly possible that the essential features of RSB might survive even in rf < 6, 
and that the appropriate analytic approach could allow to make that clear. Another 
possibility is that for c/ < 6 the droplet picture of spin glasses [33l |3l| would apply 
instead, as advocated in Ref . [5] , so that the nature of the spin glass state would change 
from being RSB like for d > 6 to being replica symmetric for d < 6. According to this 
picture Pj{q) should become just two delta functions at ±qEA in the thermodynamic 
limit, corresponding to just a state and its time reverse. (In fact, for finite systems in 
three dimensions, Pj{q) appears strikingly similar to that of the finite A^ SK mo del [35]. 
Unfortunately no systematic study has been made as to how the number of peaks, humps 
and shoulders evolves with system size; such an investigation could be very informative 
as regards the true nature of the three-dimensional spin glass state.) 

We have argued here that for the finite A^ SK model the replica symmetry breaking 
is stabilized at a finite value of K by self-energy effects. The replica symmetric state of 
the droplet picture corresponds to having K = 1. Thus were the droplet picture to be 
the valid description of spin glasses below six dimensions (and some of the authors of 
this paper would argue that this is unlikely!), the same mechanism could stabilize the 
replica symmetric state. While perturbatively there seems to be no way that the replica 
symmetric state could be stable, it is possible that if the full self-energy corrections 
about that state could be included into the calculation, then replica symmetry might 
be maintained [5]. 
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Appendix A. The A^ dependence of the self-energy 

In this Appendix we shall estimate the A^ dependence of the self-energy Hji. This 
is needed for our argument that the number of features R scale as N^^^. A direct 
calculation of Sr would be impractical: it would involve summing diagrams to all 
orders. Because of that we will obtain the A^ dependence of Hr indirectly via a study 
of the TAP equations [36] of the model. 

The TAP equations provide a non-replica way of finding single-valley correlations. 
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For the magnetization nii at site i within a single state they give 

mi = tanh /3 ^ Jijnij - /3mi ^ jfj{l - m^) + PhA . (A.l) 

The spin glass susceptibility is defined as 

.so4E(|^)^ (A.2, 

When the right-hand side is bond-averaged over the exchange interactions Jij we obtain 
Gr{{})). In the following we will only consider the case of zero magnetic field, hi = 0. 

It is convenient to express xsg iii terms of the eigenvalues of the Hessian matrix of 
the second derivatives of the TAP free energy[ 



Aij = d'^{(3FTAp)/dmidmj = -2(3'^jfjmimj - (3 Jij 

+ f/5'E4.((i-^D + (i-^n-^l%- (A-3) 



In terms of the eigenvalues A of A, 



Xsg = ^Et2' (A-4) 



N^X^ 



which in terms of the density of states p(A) becomes 



oo 

2 



XsG= lim / dXp{X)/X'. (A.5) 

The first term on the right hand side of Eq. flA.Sp is of order 1/A^ and is smaller than the 
other terms which are either of order l/N^^"^, or on the diagonal, of order 1. It will be 
dropped. (We are focusing in this work on the low-lying TAP states - the pure states 
- where the mechanism for splitting off an isolated eigenvalue as in Ref. [38] cannot 
operate.) A stable solution of the TAP equations which corresponds to a minimum 
requires all the eigenvalues of the matrix A to be positive. For pure states, p{X) is 
non-zero right to the origin; at small A (see Ref. [37]) one has that 

With this form for p(A), the integral in Eq. (lA.SP would be divergent without its lower 
cutoff at Xmin- The A^ dependence of Xmin itself can be estimated by setting 

1 = N f "^^" dX p{X) , (A.7) 

Jo 

which means that Xmin ~ N""^^^. Using this result, we can estimate the A^ dependence 
of xsG as N^^^ using Eq. (lA.Sp . Notice that this result would also apply at Tc. 

In fact there is a very simple direct argument for the behaviour at T^. According 
to Refs. [m E] for T > T^, 

XSG = ^/(^l^n , (A.8) 
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to grow faster than the predicted N^^^. This is due to finite size effects, see text. 
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Figure A2. Scaling plot of the averaged integrated density of states of the Hessians 
for different system sizes. 



SO that as |t| goes to zero, that is, at Tc, xsg ~ N^^^. 

We would not expect that bond-averaging xsg to get Gij(O) will modify this A^ 
dependence since single-valley quantities are expected to be self-averaging. 

Thus for T < Tc, the single valley spin glass susceptibility Gr{0) diverges as N^^^. 
This implies that the typical value of K, the order of replica symmetry breaking in a 
finite system of A^ spins, will be via Eq. ([6]) of order tN^^^. The data in Fig. [T]is clearly 
consistent with this expectation. 
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We have tested numerically the prediction xsg ~ N^^^ using the iteration procedure 
described in Ref. [39]. This method allows us to find many TAP states even for large 
system sizes (see Ref. [39] for a discussion of the proximity of the iteration algorithm to 
a dynamical critical point and the infiuence of this on the free energy of the states found 
- here we chose the proximity for each system size in such a way that we get the same 
free energy range for all system sizes). Using systems with N = 100, 200, 283, 400, 566 
and 800 at a temperature of T = 1.0//3 = 0.2, we have calculated the Hessian for every 
state found. The diagonalization of the Hessians was done using arbitrary precision 
arithmetic (with an accuracy of 50 decimal digits) since these matrices are extremely 
ill-conditioned and so they are hard to diagonalize: standard packages, for instance the 
LAPACK routines, fail at the task. The eigenvalues were used to calculate xsg as in 
Eq. ( ]A.4[) . The results are shown in Fig. lAll 

Surprisingly, xsG appears to grow faster than N^^^. We believe, however, that this 
is a finite size effect. To back up this claim, we show in Fig. IA2I a scaling plot of the 
averaged integrated eigenvalue density -D(A) of the Hessians of sizes N = 200,400 and 
566. 

The expectation is that in the thermodynamic limit this function goes as -D(A) = 
-Doo(A) ~ A'^/^ for small A (corresponding to p(A) ~ A^^^). For finite TV there is a 
cutoff around A ~ N""^^^. The natural expectation is that this cutoff is of the form 
D{X) = Doo(A)/(iVA3/2) ^i^i^ a scaling function f{x). This is verified in Fig. ESI The 
arguments sketched above which lead to the prediction xsg ~ A^^^^ for large A^ can 
only be expected to be valid when the interval in which -D(A) ~ A'^^^ prevails is large 
enough. This interval can be identified as the horizontal (or nearly horizontal) stretch 
in Fig. IA2[ Clearly, this interval is very small as it is not even one decade for the 
available system sizes. The conclusion is, therefore, that we cannot yet expect to see 
the asymptotic scaling behaviour. It is not possible to go to larger system sizes as 
the arbitrary precision diagonalization of the Hessians becomes computationally too 
expensive. 

Appendix B. Free energy fluctuations and J{p) 

According to Ref. [10], the exact expression for the quantity J{p) for a finite number of 
replica symmetry breaking steps K and for the truncated model with Hamiltonian 

^ = "2 E ?a/3 - -^ E 9"/39/37?7« - Y^ E ?a/3 (B-l) 

a,/3 cri/3,7 a, (3 



IS 



K+1 

J{p) = -Y^ /io(A;)/io(0 log(p' + A(0; k, I)) , (B.2) 

k,l=l 



where 



1 1 



k > r — 1 



/..(fc) ={ -^-- ^^ , (B.3) 

Pr + l 
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Pk = %K^ , (B.4) 

w K 

A(r; fc, /) = 2y^ (i(fc - 1)^ + ^(/ - 1)^ - ^' - ^) ' (B-5) 

and qK is the solution of 



t - wqk + y[l- -^j qK = 0- (B.6) 

When the expressions for fir{k) and /ir(/) are inserted, Eq. (1B.2P can be rewritten as 

,%PkPi''^[p^ + x{o-k + i,i) p2 + x(o-kj + i) )■ ^-^ 

2 

We are interested in the behaviour for small p where J(p) diverges as p approaches |^. 
It is easy to see that the first two terms of the former expression are well behaved in 
this limit. The divergence in p must therefore come from the last term, which will be 
denoted by J{p). Defining 



„2 _ KV 1 

yql 3 



a;2 = —^ - - (B.8) 



and renumbering the sums to start from it can be cast in the form 



wHx^+^^ i^-l 



j^p)=-^j-r^i: 



^yp^ ,% (fc + I) (/ + 1 



X log ( -^ + ^^ + ^^ X^ + (fe+l)- + (/ + in (B g. 

While p was a finite-dimensional wave vector in Ref . [10] , here we consider it as a proxy 
for the self-energy as we are dealing with the SK model. Substituting Sr = cN~^^^ 
for p"^ and making the usual replacement K = c'tN^'^ (where the constant c' is large 

enough to guarantee stability) yields (to leading order in A^) x^ = ^-^ | and 

J = SF'^ = N^^'^fit). The function f{t) is defined by the remaining prefactors and 
the sums (with upper bounds set to infinity) in Eq. ( IB. 91) . This shows that the typical 
sample-to-sample fluctuations of the free energy are of order A^^/®. 
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